function y_psd = lombpsd(x_time, y_data, x_freq);

% LOMBPSD - use Lomb's method to compute normalized PSDs
%
%    PSD = LOMBPSD(T, X, F) uses Lomb's method to compute normalized PSDs
%    evaluated at frequencies F from the vector X sampled at times T to
%    obtain statistical measures of PSD features (Numerical Recipes,
%    section 13.8).
%
%    Note that this routine only works for individual data sets.  If X and T
%    are matrices, they are reshaped into vectors and one PSD is calculated.

% By:   S.C. Molitor (smolitor@eng.utoledo.edu)
% Date: January 4, 2001
%
%    Note that there are two implementations of this code in LOMBPSD.M.  The
%    first is straight from equations 13.8.4 and 13.8.5.  The second uses
%    trigonometric identities to reduce the number of calls to SIN() and
%    COS(), which are slow.

% initialize output

y_psd = 0;

% validate arguments

if (nargin ~= 3)
   msgbox('Invalid number of arguments', 'LOMBPSD Error', 'warn');
   return
elseif (~isnumeric(x_time) | isempty(x_time))
   msgbox('T must be a numeric vector', 'LOMBPSD Error', 'warn');
   return
elseif (~isnumeric(y_data) | isempty(y_data))
   msgbox('X must be a numeric vector', 'LOMBPSD Error', 'warn');
   return
elseif (~isnumeric(x_freq) | isempty(x_freq))
   msgbox('F must be a numeric vector', 'LOMBPSD Error', 'warn');
   return
elseif (prod(size(x_time)) ~= prod(size(y_data)))
   msgbox('T and X must have the same number of elements', 'LOMBPSD Error', 'warn');
   return
end

% convert data & time into zero-mean 1 x N row vectors

t = reshape(x_time, 1, prod(size(x_time)));
tm = mean(t);
t = t - tm;
x = reshape(y_data, 1, prod(size(y_data))) - mean(mean(y_data));
xm = mean(x);
x = x - xm;

% convert frequency to M x 1 column vector (rad/sec)
% convert data to M x N array (identical rows)
% create M x N array of frequency * time

w = 2 * pi * reshape(x_freq, prod(size(x_freq)), 1);
h = ones(size(w)) * x;
wt = w * t;

% Method 1: straightforward evaluation of eqns 13.8.4 & 13.8.5
% has two calls to SIN() and COS() on M x N arrays

% calculate M x 1 product of angular frequency and tau
% create M x N array of frequency * (time - tau)
% calculate sine & cosine of this M x N array
% create M x N array of data * sine or cosine
% sum across columns to create M x 1 PSD vector
% PSD is transposed back to row vector (1 x M)

%wtau = 0.5 * atan2(sum(sin(2 * wt), 2), sum(cos(2 * wt), 2));
%arg = wt - (wtau * ones(size(t)));
%coswt = cos(arg);
%sinwt = sin(arg);
%hcoswt = h .* coswt;
%hsinwt = h .* sinwt;
%y_psd = ((sum(hcoswt, 2) .^ 2) ./ sum(coswt .^ 2, 2) + (sum(hsinwt, 2) .^ 2) ./ sum(sinwt .^ 2, 2))';


% Method 2: evaluate eqns 13.8.4 & 13.8.5 with trig identities
% has single call to SIN() & COS() on M x N arrays

% calculate sin(wt) & cos(wt)
% calculate product of data and sin(wt) & cos(wt)
% calculate products & squares of sin(wt) & cos(wt)

coswt = cos(wt);
sinwt = sin(wt);
hcoswt = h .* coswt;
hsinwt = h .* sinwt;
cos2wt = coswt .^ 2;
sin2wt = sinwt .^ 2;
cossinwt = coswt .* sinwt;

% sum arrays across columns (creates M x 1 vectors)

shcoswt = sum(hcoswt, 2);
shsinwt = sum(hsinwt, 2);
scos2wt = sum(cos2wt, 2);
ssin2wt = sum(sin2wt, 2);
scossinwt = sum(cossinwt, 2);

% calculate PSD from M x 1 sum vectors
% PSD is transposed back to row vector (1 x M)

wtau = 0.5 * atan2(2 * scossinwt, scos2wt - ssin2wt);
coswtau = cos(wtau);
sinwtau = sin(wtau);
shcos = coswtau .* shcoswt + sinwtau .* shsinwt;
shsin = coswtau .* shsinwt - sinwtau .* shcoswt;
s2cos = (coswtau .^ 2) .* scos2wt + 2 * sinwtau .* coswtau .* scossinwt + (sinwtau .^ 2) .* ssin2wt;
s2sin = (coswtau .^ 2) .* ssin2wt - 2 * sinwtau .* coswtau .* scossinwt + (sinwtau .^ 2) .* scos2wt;
%y_psd = (((shcos .^ 2) ./ s2cos) + ((shsin .^ 2) ./ s2sin))';

% prevent division by zero in PSD calculation
% substitute mean value where either denominator is zero
% NOTE:  this is a departure from Numerical Recipes algorithm!

z2cos = (s2cos == 0);
z2sin = (s2sin == 0);
%z2cs = (z2cos | z2sin);
r2cos = (1 - z2cos) .* (shcos .^ 2) ./ (s2cos + z2cos);
r2sin = (1 - z2sin) .* (shsin .^ 2) ./ (s2sin + z2sin);
%y_psd = ((r2cos + r2sin) .* (1 - z2cs) + xm * z2cs)';
y_psd = (r2cos + r2sin)';
return

